active_sites = [1 2 4 7 8];
passive_sites = setdiff(1:9,active_sites);
control_sites = 10;

figure;hold on
mean_resp_by_site = [];

for cur_exp = 1:9
    rps_resp_cnt = 0;
    rps_resp = [];
    for time_point = 1:4
        for layer = 1:4
            cur_act = proj_meta(cur_exp).rd(layer,time_point).act(:,1:7500);

            % generate binary activity during feedback.
            velM_bi = proj_meta(cur_exp).rd(layer,time_point).velM_smoothed(1:7500);
            velM_bi = -velM_bi; % make the velocity of forward-running positive.
            velM_bi(velM_bi > 0.001) = 1;
            velM_bi(velM_bi <= 0.001) = 0;
            
            % generate the frames of ps onset.
            ps_ons = find(diff(proj_meta(cur_exp).rd(layer,time_point).ps_id(1:7500))>1);
            % remove early & late ps.
            ps_ons(ps_ons > 7425) = [];
            ps_ons(ps_ons < 6) = [];
            
            rps_nbr = 0; % number of real ps
            rps_id = [];
            for psid = 1:length(ps_ons)
                if sum(velM_bi(ps_ons(psid)-5:ps_ons(psid)+15)) == 21
                    rps_nbr = rps_nbr+1;
                    rps_id(rps_nbr) = psid;
                end
            end
            
            for ROI_id = 1:size(cur_act,1) %mm_cells{cur_exp,layer}
                for rps_act_nbr = 1:length(rps_id)
                    rps_frame_id = ps_ons(rps_id(rps_act_nbr));
                    rps_resp_cnt = rps_resp_cnt+1;
                    rps_resp(rps_resp_cnt,:) = cur_act(ROI_id,rps_frame_id-7:rps_frame_id+75);
                    resp_by_ps(cur_exp,time_point,rps_act_nbr) = mean(rps_resp(rps_resp_cnt,23:30))-mean(rps_resp(rps_resp_cnt,1:7));
                end
            end
        end
    end
    mean_resp = mean(rps_resp);
    mean_resp = mean_resp - mean(mean_resp(1:7));
    mean_resp_by_site(cur_exp,:) = mean_resp;
    if sum(cur_exp == active_sites)
        plot(mean_resp)
    else
        plot(mean_resp,'r')
    end
end

sem_resp_act = std(mean_resp_by_site(active_sites,:))/sqrt(size(mean_resp_by_site(active_sites,:),1));
sem_resp_pas = std(mean_resp_by_site(passive_sites,:))/sqrt(size(mean_resp_by_site(passive_sites,:),1));

figure;hold on
area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean(mean_resp_by_site(active_sites,:))+sem_resp_act fliplr(mean(mean_resp_by_site(active_sites,:))-sem_resp_act)])
area([-7/15:1/15:75/15 75/15:-1/15:-7/15],[mean(mean_resp_by_site(passive_sites,:))+sem_resp_pas fliplr(mean(mean_resp_by_site(passive_sites,:))-sem_resp_pas)])
plot([-7/15:1/15:75/15],mean(mean_resp_by_site(active_sites,:)))
plot([-7/15:1/15:75/15],mean(mean_resp_by_site(passive_sites,:)),'r')
xlim([-0.5 5])

figure;hold on
plot(smooth2(squeeze(mean(mean(resp_by_ps(active_sites,:,:),1),2)),3))
plot(smooth2(squeeze(mean(mean(resp_by_ps(passive_sites,:,:),1),2)),3),'r')